LAD calibration
Nathéo BEAUCHAMP
2025-09-24
Motivation of the study
This study aims to calibrate Leaf Area Density (LAD), which is defined as the area of leaves per unit of tree crown volume (in \(m^2.m^{-3}\)). This parameter forms the basis of numerous process-based forest models, such as physiological models that consider photosynthesis and transpiration (e.g. SurEau by Cochard et al., 2021), and forest gap models that estimate light interception (e.g. ForCEEPS by Morin et al., 2021, and Samsara2 by Courbaud et al., 2015).
Although the leaf area density can be directly estimated from a sample of trees (Stadt & Lieffers, 2000), this method is laborious and difficult to replicate for every tree in a stand. Therefore, indirect methods can be employed, such as model inversion using a hemispherical photograph (Courbaud et al., 2003) or the more recent LiDAR methodology (Wei et al., 2020). However, there is a significant lack of research on LAD in the scientific literature, as there are no large-scale databases containing values for a wide range of European species. Furthermore, the few reported species-specific values vary considerably between studies (Ligot et al., 2014).
The aim of this study is to estimate leaf area density (LAD) for nine different species/groups: Fagus sylvatica, Picea abies, Quercus sp., Carpinus betulus, Pseudotsuga menziesii, Pinus sylvestris, Abies alba, Larix decidua and ‘other species’. Additionally, Nock et al. (2008) demonstrated intraspecific variation in LAD, with intracrown leaf area decreasing by up to 40% with increasing tree age. This led us to consider the effect of diameter as well. Finally, as crown dimensions are very sensitive to local competition (Touzot et al., 2025), we hypothesise that LAD may also vary in response to competition.
Methodology
In this study, we used a model inversion methodology, using the ray-tracing model SamsaraLight (Courbaud et al. 2003). We want to estimate the SamsaraLight tree-level parameter LAD using measurements of transmitted light \(PACL_{obs}\) (for percentage of above light canopy), based on hemispherical photographs, within 45 plots of different composition and structures. For each plot, trees are inventoried with their crown dimensions, allowing to derive an estimated transmitted light \(PACL_{est}\) on virtual sensors with the ray-tracing model SamsaraLight. By doing so, for each virtual sensor of each plot, we have the observed transmitted light \(PACL_{obs}\), and the estimated transmitted light given a tree-level value of LAD \(PACL_{est}\).
The aim of the calibration is to find a tree-level LAD model that minimize the error between \(PACL_{obs}\) and \(PACL_{est}\). To do so, we defined the LAD for a tree \(i\), of species \(s\) in a plot \(p\) as \(LAD_{i,s,p} = a_{0,s} + a_{0,p} + a_{1,s}.DBH_i + a_{2,s}.BAH_i + \epsilon\), where \(a_{0,s}\) is a species-specific intercept, \(a_{0,p}\) is a random hierarchical site effect, nested within an origin effect (4 origins of the plots within UCLouvain, INRAe, Cloture and IRRES), \(DBH_i\) is the diameter at breast height of tree \(i\) , with the associated estimated species-specific coefficient \(a_{1,s}\), \(BAH_i\) the competition index of the tree \(i\) (i.e. the sum of basal area of trees higher than the tree \(i\)), with the associated estimated species-specific coefficient \(a_{2,s}\), and \(\epsilon\) the normal error term. We used a Bayesian approach to estimate the posterior distribution of parameters that minimize the log-likelihood between observed and simulated PACL. As computation of \(PACL_{est}\) needs an external call of SamsaraLight model, we used the R package ‘BayesianTools’ to run a MCMC sampling of parameters with the “DREAMzs” sampling algorithm.
For this study, we used the SamsaraLight model with the R package SamsaraLight (Beauchamp et al., work in progress, intended for publication in the Journal of Open Source Software). The SamsaraLight model was initially developed in Java and included in the Capsis simulation platform (SamsaraLightLoader by Gauthier Ligot). It could be called from R software via the “RCapsis” R package (Fortin et al., 2021). To speed up calculation times and make SamsaraLight easier to use in an R environment, we created the “SamsaraLight” R package. The source code is written in C++ to enable fast calculations and communication with the R environment is handled by the Rcpp package (Eddelbuettel & Balamuta, 2018). The “SamsaRaLight” R package reduces calculation time compared to the previous method (the RCapsis R package calling the Java model SamsaRaLightLoader) as it does not require the creation of a Java server nor Java objects. It also facilitates the use of SamsaraLight within R, as it does not require the creation of an inventory file, instead using R objects and functions directly. The package is functional and stored in a GitHub repository, but it is still subject to private restrictions as discussions regarding co-authorship and the intellectual property rights of the source code need to be concluded.
Calibration plots
The 45 calibration plots comes from 4 different sources we wiil call IRRES (9 plots, Gauthier Ligot), CLOTURE (23 plots, Gauthier Ligot), INRAe (10 plots, Philippe Balandier) and UCLouvain (3 plots, Mathieu Jonard).
IRRES
CLOTURE
INRAe
UCLouvain
Preliminary analysis
As a preliminary analysis, we tried to fit a mean LAD value for each sensor of each site. To do so, we run the SamsaraLight model on each stand by fixing the same crown LAD values of each tree, and for each sensor, we computed the residuals between the estimated PACL value on the virtual sensor, and the observed PACL on the field for this sensor. We did so for 500 values of LAD from 0.01 to 5. By doing so, we could find for each sensor the LAD that minimizes the residuals.
Many sensors did not converged (i.e. best fitting LAD was 5 as the residuals function showed an asymptotic form). Indeed, we could increase as much as we want the trees’s LAD, the estimated PACL could not decrease as much as the observed one. PACL is principally linked to the crown geometries, and secondary by the crown LAD. Thus, there are two cases where crown LAD did not influence PACL: (i) for open stands, where the high light on the ground comes principally from unobstructed rays, (ii) for highly crowded stands, where the very low light on the ground comes from the high number of intercepted crowns. Another reason could be linked to the wrong representation of leaves and crown above the sensor, where a simple branch can obstruct the sensor in the field where we considered a simple crown shape that result in incapabilities of representing such a light obstruction with the SamsaraLight model. Consequently, we remove the sensors that did not converged (465/1121 sensors), leading to remove consequently 3 sites where all sensors did not converged: Cloture11, Cloture15, Cloture2.
The mean LAD over all sensors that converged was about 0.944
Bayesian calibration
Clean model: species X DBH with covariance
100,000 iterations and 3 chains
## id_model n_iterations n_analysis filter_sensors site_rd_effect
## 1 1 2e+05 5000 TRUE TRUE
## 2 2 2e+05 5000 TRUE TRUE
## 3 3 2e+05 5000 TRUE TRUE
## origin_rd_effect intercept_per_sp intercept_sp_pooling dbh_effect dbh_per_sp
## 1 TRUE FALSE FALSE TRUE FALSE
## 2 TRUE TRUE TRUE TRUE TRUE
## 3 TRUE TRUE TRUE TRUE TRUE
## dbh_sp_pooling consider_covariance compet_effect
## 1 FALSE FALSE FALSE
## 2 TRUE FALSE FALSE
## 3 TRUE TRUE FALSEComputation time
## model chain time_hours
## 1 1 1 310.8808
## 2 1 2 317.5449
## 3 1 3 300.6420
## 4 2 1 319.6946
## 5 2 2 317.3686
## 6 2 3 306.4468
## 7 3 1 317.5156
## 8 3 2 301.6323
## 9 3 3 317.2304Observe convergence
Effect of DBH for each species
## # A tibble: 27 × 6
## # Groups: model, chain [9]
## model chain subchain rho_lower rho_median rho_upper
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 1 1 1 0 0 0
## 2 1 1 2 0 0 0
## 3 1 1 3 0 0 0
## 4 1 2 1 0 0 0
## 5 1 2 2 0 0 0
## 6 1 2 3 0 0 0
## 7 1 3 1 0 0 0
## 8 1 3 2 0 0 0
## 9 1 3 3 0 0 0
## 10 2 1 1 0 0 0
## # ℹ 17 more rows
Effect of BAxDBH
100,000 iterations and 1 chain
## id_model n_iterations n_analysis filter_sensors site_rd_effect
## 1 1 1e+05 5000 TRUE TRUE
## origin_rd_effect batot_in_site_effect intercept_per_sp intercept_sp_pooling
## 1 TRUE TRUE FALSE FALSE
## dbh_effect dbh_interaction_batot dbh_per_sp dbh_sp_pooling
## 1 TRUE TRUE FALSE FALSE
## consider_covariance compet_effect
## 1 FALSE FALSEComputation time
## id_model time_hours
## 1 1 54.08484Observe convergence
Effect of DBH X batot on a mean species
Output analysis
Observe initial data
##
## Call:
## lm(formula = batot_m2ha ~ order, data = data_output_tree)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.733 -7.807 0.806 5.892 34.167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.4670 0.1306 164.38 <2e-16 ***
## ordergymnosperm 11.6324 0.1889 61.59 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.527 on 8167 degrees of freedom
## Multiple R-squared: 0.3172, Adjusted R-squared: 0.3171
## F-statistic: 3794 on 1 and 8167 DF, p-value: < 2.2e-16Stand variables
Crown dimensions
Tree-level consequences on tree leaf area
Stand-level consequences on leaf area index LAI
Consequences on light
Tree-level
Stand-level
